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Abstract 

We formulate a stochastic equation to model the erosion of a surface with fixed incli- 
nation. Because the inclination imposes a preferred direction for material transport, the 
problem is intrinsically anisotropic. At zeroth order, the anisotropy manifests itself in a lin- 
ear equation that predicts that the prefactor of the surface height-height correlations depends 
on direction. The first higher-order nonlinear contribution from the anisotropy is studied 
by applying the dynamic renormalization group. Assuming an inhomogeneous distribution 
of soil substrate that is modeled by a source of static noise, we estimate the scaling expo- 
nents at first order in e-expansion. These exponents also depend on direction. We compare 
these predictions with empirical measurements made from real landscapes and find good 
agreement. We propose that our anisotropic theory applies principally to small scales and 
that a previously proposed isotropic theory applies principally to larger scales. Lastly, by 
considering our model as a transport equation for a driven diffusive system, we construct 
scaling arguments for the size distribution of erosion "events" or "avalanches." We derive a 
relationship between the exponents characterizing the surface anisotropy and the avalanche 
size distribution, and indicate how this result may be used to interpret previous findings of 
power-law size distributions in real submarine avalanches. 
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1 Introduction 



One of the great challenges faced by modern research in nonlinear physics is the construction 
of predictive theories for systems in which the underlying equations of motion are not known. 
An example that has recently created much interest is the flow of sand, or, more generally, 
granular fluids ||], |2|. One expects that the Navier-Stokes equations do not apply to granular 
flow because of the dissipative nature of grain-grain collisions. Although physical arguments 
have been used to deduce equations of motion applicable to certain situations (e.g., see Ref. ||), 
a comprehensive theory of granular flow remains elusive. 

Still more complicated than sand is the flow of wet sand j|. Viewed naively, in such a 
situation we retain the complexity of the Navier-Stokes equations and add to it the complicated 
frictional stresses of a granular heap. In this paper we shall concern ourselves with wet sand 
in the form of a particular problem in geomorphology. Specifically, we study the erosion of a 
landscape due, usually but not exclusively, to the flow of water over it. Despite the obvious 
difficulties of this problem, our aim is to obtain some simple results that offer fundamental 
insight into erosion. Why we expect to do so deserves some further comment. 

There is now a wealth of empirical evidence that shows that real landscapes exhibit some 
form of scale invariance [||. These scaling laws come in many forms. Perhaps best known are 
those that describe the branching of river networks ||. These are not the scaling laws that 
interest us here, however. Instead, we wish to view the problem of topography more generally, 
by examining statistical properties of surfaces h(x), where x denotes the horizontal position and 
h is the topographic height or elevation. 

The statistics that primarily interest here are the height-height correlation functions 

C(r) = (|Mx + r)-Mx)| 2 )i/ 2 . (I) 

Scale invariance comes in the form of self-affinity ||. In other words, C(r) ~ r a , where a is 
known as the roughness exponent. Various empirical measurements in different sorts of terrain 
show that this power-law form may hold over an order of magnitude or more, with some mea- 
surements indicating that a is small (0.30 < a < 0.55) [7-13] while others show it to be large 
(0.70 < a < 0.85) [9-16]. 

Although there is not much agreement on the value of the scaling exponent a, the occurrence 
of scaling itself is fairly common. We are therefore led to consider theoretical models whose 
solutions also exhibit scaling. These models are stochastic partial differential equations that are 
Langevin equations for the evolution of a surface @ . One of our goals here is to derive such an 
equation that predicts some aspects of the observed scaling. Our hope is that our predictions 
are general and independent of details such as material properties, climate, etc. Thus we hope 
that our model exhibits some degree of universality 

For our model to exhibit universality, we must identify a class of topographic evolution 
problems for which we may make quantitative predictions. The class of problems we discuss 
here are problems in which symmetry is broken by the existence of a preferred direction — 
downhill — for the flux of eroded material. Following work we have already reported in a brief 
Letter M, we derive an anisotropic noisy diffusion equation to describe erosion at the small 
length scales where the preferred direction is fixed throughout space. The linear regime of this 
equation predicts that C(r) is anisotropic at the level of its prefactors. The predicted anisotropy 
is testable, and empirical studies in progress show that it works with unusual generality |ict |. 
Under the additional assumptions that the flux of eroded material increases with increasing 
distance downslope and that the dominant effects of noise are fixed in space, we find, using 
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the dynamic renormalization group, that not only is C(r) anisotropic, but that it scales with 
different exponents that correspond to the downhill direction and the direction perpendicular 
to the downhill direction. This result is also testable, and we present one example, made from 
the topography of the continental slope off the coast of Oregon, in good agreement with our 
predictions. 

An additional conclusion of our study concerns the wide range of values of a that have been 
reported in the literature. It has been proposed previously [11| that observations of a ~ 0.4 



could be explained by the Kardar-Parisi-Zhang (KPZ) equation [12 



^ = ^ h+ ^ Vh f + v . (2) 

In its application to geomorphology, u is a topographic diffusivity coefficient, A is related to 
the velocity of erosion in the direction normal to the surface, and r\ is a source of random 
noise that is uncorrelated in space and time. We observe here that if equation (|2|) really does 
capture some aspects of topographic evolution, then it applies only to those cases in which a 
is found to be small. It turns out that most observations of small a are made at large length 
scales where no preferred direction is easily identified [7-13], whereas observations of large a 
are usually associated with small length scales [9-16]. Because the average results predicted 
by our anisotropic theory are consistent with these large-a observations, we can tentatively 
identify two "universality classes" of topographic evolution. In the KPZ class, topography 
evolves isotropically (perhaps due to internal tectonic stresses) at large length scales and yields 
small roughness exponents, while in our anisotropic class, topography evolves erosively, "one 
slope at a time," at small length scales and yields large roughness exponents. We lend some 
support to this conclusion by showing evidence of such a crossover in a single topographic 
dataset. 

The general framework of our theory also allows us to make some contact with the larger 
field of "self-organized criticality" (SOC) [jl3| . Specifically, sloping submarine topography gives 
rise to underwater avalanches. These avalanches create flows, known as turbidity currents fl4| ], 
that eventually come to rest as sedimentary deposits called turbidites [15]. A number of recent 



studies have indicated that the size distribution of these natural avalanches may follow the 



power-law scaling predicted by SOC sandpile models |16], [17], 18]. Here we show how our theory 
for topographic evolution may be linked to the SOC theory for avalanche sizes. Specifically, we 
derive a relation between the anisotropic correlations of the slope and the size distribution of 
the avalanches. 

This long introduction will have succeeded if the reader is convinced that the concepts of 
scaling and universality may have some applicability to understanding some generic features of 
our natural environment. In this spirit, it is our pleasure to dedicate this paper to Leo Kadanoff. 



2 Mathematical formulation 

We begin our discussion with a brief introduction to stochastic equations for surface growth 



19|, |20|]. We first review some standard isotropic models, and briefly remark on their applicability 
to geomorphology. We then introduce our anisotropic model. Analysis of the model is deferred 
to the following section. 
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2.1 Isotropic surface growth 



Our objective is to determine the evolution of the surface h(x,t), where, as we have already 
stated, h is the height of the surface at position x and time t. We assume that h is single 
valued — that is, overhangs are not allowed. The general form of an equation for h that we 
consider here is 

^M=^(x,*)]+^x,*). (3) 

J- represents the flux of eroded material, and r\ is a source of random noise that allows us to 
include random fluctuations in the growth process. In the absence of specific information on T, 
one generally seeks to first identify all applicable physical symmetries and conservation laws. 
This then allows the construction of the simplest possible form of T compatible with these 
constraints [pCj ]. 

Of the vast number of equations such as (^) that have been proposed in recent years [2^] , 
here we restrict the discussion to models that have been proposed for the study of erosion at 
large length and long time scales. These may be divided roughly into two categories: those which 
conserve a material flux J and those which do not. In the conservative models, T = — V • J, 
where J is the current of material. The simplest of these has J = — V/i and no noise, leading to 
the classical diffusion equation, which, in geomorphology, was first popularized by Culling |)21|: 



The diffusion equation alone, however, cannot explain the formation of stable self-affine land- 
scapes. If we add uncorrelated noise, we obtain the so-called Edwards- Wilkinson equation 



[22, With this equation we can obtain true self-affine surfaces, but in the relevant number 



of dimensions (i.e., when the dimensionality of the position vector x is d = 2), the noisy version 



of equation (Q) predicts that correlations decay logarithmically (i.e., a = 0) |22|, 23 1. Thus nei- 
ther the deterministic nor stochastic form of equation is compatible with the aforementioned 
observational evidence. 

Partly as a remedy for this problem, non-conservative equations have been proposed. The 
most popular of these is the KPZ equation (0). As Sornette and Zhang |11] have remarked, 



the KPZ equation is attractive as a model of erosion because it is the simplest surface growth 
equation that can generate a nontrivial (a ^ 0) self-affine landscape. Specifically, the KPZ 
equation contains the necessary ingredients of nonlinearity and noise. Roughly speaking, the 
nonlinearity results from a uniform rate of erosion, at all locations x, in the direction normal to 
the surface, and the noise accounts for the irregularities of the process in time and space. The 
roughness exponent reported for KPZ in d = 2 by numerical integration varies between 0.2 and 



0.4 [24, 25, 26]. As discussed above, these values are in reasonable agreement with observations 
at large length scales, where one finds that 0.30 < a < 0.55 [7-13]. Thus the KPZ scenario of 
non-conservative isotropic growth normal to the surface may indeed apply to some aspects of 
large-scale landscape evolution. 



2.2 An anisotropic model 

As discussed in the introduction, the predictions of the KPZ equation do not agree with many 
measurements made on landscapes at small length scales [9-16]. Thus some other physical 
mechanisms must be dominant in this range. Here we propose that evolution at small length 
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Figure 1: Schematic configuration of an anisotropic landscape in d = 2. 

scales is strongly influenced by the breaking of symmetry induced by the presence of a slope of 
fixed inclination. 

Figure [l] depicts the framework for our theory. The vector is the growth direction in our 
parametrization. Note that is measured from the top of the slope downwards. The action 
of gravity selects a preferred direction given by the vector en, which is essentially defined by 
the direction "downhill," and which coincides with the average direction of the flow of material. 
The symbol e_|_ stands in general for the subspace of directions perpendicular to en. For a real 
landscape, this subspace is spanned by a single vector. Later we generalize to landscapes on 
]R rf , in which case the perpendicular subspace has dimension d — 1. The framework of Fig. [l| is 
completed by imposing fixed boundary conditions at the top of the slope (i.e., at x\\ = 0), or by 
imposing the symmetry xii — > —%\\- 

Because Fig. [l] explicitly distinguishes between the two directions ej_ and en, we expect this 
same anisotropy to be reflected in the correlation function C(r). Thus, if h is self-affine, we 
expect different roughness exponents for correlations measured in each of the directions en and 
e_|_. Specifically, we define an and a± such that C||(x||) ~ |^||| a|1 for correlations along fixed 
transects Xj_ = const., and Cj_(xj_) ~ |xj_| a± for correlations along fixed transects xjj = const., 
where in general aii / a±. These relations are summarized by the single scaling form 

C(x||,x ± ) ~ & a IIC(& -1 S||,6- c «xi), (5) 

where £y is the anisotropy exponent. The roughness exponents ay and a± are related to £|| by 

an 

a± = — . (6) 

Ml 

The anisotropy exponent <ji accounts for the different rescaling factors along the two main 
directions. Since the space is anisotropic, when performing a scale change, we must rescale xu 
and x^ by different factors bn and b±, respectively, if we are to recover a surface with the same 
statistical properties. We assume that this scaling is self-affine, such that 

C|| = "j — r- = const. (7) 

logon 

Note that 0i defines only the ratio of the roughness exponents, but not their precise magnitudes; 
moreover, the scaling form (|5|) is not unique. We can also express the rescaling along the direction 



4 



ej_ by writing 

C(s||,xj_) ~ b a± CQ>-t x x\\,b- 1 x±), (8) 

where we have used the anisotropy exponent £j_ = 1/C||- Both scaling forms (||) and (§) are 
completely equivalent. 

We seek a single stochastic equation for the landscape height h. We assume that the deter- 
ministic motion of the underlying soil is locally conserved such that 

f = -V. J+ ,, (9) 

where J is the current of soil per unit length. The soil however is not globally conserved, since 
it is lost at the bottom boundary. Conservation is also broken by the addition of a stochastic 
noise term 77, discussed below. 

Physically, the current J is expected to reflect two effects. First, we expect a local isotropic 
diffusing component, tending to smooth out the surface. Second, we expect an average global 
flow of dragged soil, directed mainly downhill. Thus we postulate the following form for the 
current: 

J = —vVh - YV\\h. (10) 

The first term corresponds to Fick's law for diffusion, and represents the isotropic relaxational 
dynamics of the soil. The second term represents the average flow of soil that is dragged downhill, 
either due to the flow of water or the scouring of the surface by the flow of the soil itself. The 
direction of this term is given by the vector Vn/i = d»h en. The term V plays the role of an 
anomalous anisotropic diffusivity. In order to gain insight into the role of T, consider the case 
in which erosion results from the stress exerted on the soil bed by an overland flow Q of water, 
where Q is the volumetric flow rate though unit area perpendicular to the direction of steepest 
descent. The greater Q is, the stronger is the stress ||, ^7|. Moreover, since Q flows downhill, it 
increases with distance downslope. Thus T must be an increasing function of xn. Since the fixed 
inclination implies that h increases with xn, we choose to parameterize the anomalous diffusion 
as a function of the height such that T = T(h). Defining T(h) = Xq + g(h), with g(0) = and 
G(h) = J g(h)dh, we substitute Eq. flu]) into (|). Then, since 

g(h)d ]l h = ^^-d ]l h = d ll G(h), (11) 

where we have used the chain rule for the second equality, we obtain 

^t = Ul] df l h + v ± Vlh + df ] G(h)+r), (12) 

where v± = v and vn = v + Aq. 

We may advance still further by making some additional assumptions. Assuming that T(h) 
is an analytical function, we can perform a Taylor expansion in powers of h. Since all odd powers 
of h must vanish in order to the preserve the joint symmetry h — > — h, J — > —3 in Eq. @, we 
are left at lowest order with g{h) ~ A2/1 2 . By dimensional analysis (see next section and the 
Appendix) one can check that all the terms h 2k in this expansion are relevant under rescaling. 
However, the flux Q{xn) of the erosive agent (water or soil) flowing on the surface should grow 
no faster than Q(xu) ~ xjf. Then, taking h ~ xh, we find that the terms in g(h) should be of 
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order h d or less. Specializing to the case of d = 2 (i.e., real surfaces), we then find it reasonable 



to truncate g at second order. Equation (12) then takes the form 



^ = v$h + y^h + \tj{{h*)+i i , (13) 

where A = A2. 

Equation (|l^) constitutes our full nonlinear theory. Note that it differs significantly from the 
anisotropic driven diffusion equation of Hwa and Kardar [^, ^] . The differences are essentially 
due to the form of our current J, which in our case is suggested not only from symmetry 
principles, but also from the physics of erosion. We note additionally that, unlike some previous 
models of landscape erosion that couple two equations — one for the landscape h, and one for 
the overland flow Q [41-44] — here we have implicitly assumed that the effects of Q may be 
subsumed into the functional dependence of J on h. Specifically, our fundamental assumption 
of a preferred direction for the current J can be traced to global constraints imposed by the 
fixed inclination. Thus average effects of Q — for example, the fact that the erosion rate increases 
with X|| for A2 > — survive our local formulation. 

It remains to discuss the issue of noise. We distinguish two possible different sources. First, 
we may allow a term of annealed (time-dependent) noise, rj(x, t), depending on time and position, 
and describing a random external forcing due, for example, to inhomogeneous rainfall. We 
assume that this noise is isotropic, Gaussian distributed, with zero mean, and correlations 

( I) (x,t) t) (x',(')) = 2D5 W {x-x')5(t-t'). (14) 

Second, we may have a term of quenched (time- independent) noise to account for the hetero- 



geneity of the soil, mimicking the variations in the erodibility of the landscape [30]. The notion 
of quenched noise is common in the study of interface growth in a disordered medium, close to 
the depinning transition [pi]] . In this case, the noise is a function of position and height, r/(/i,x), 
with correlations given in general by 

(r]{h, x.)r]{ti, x')> = 2A(h - h')5 W (x - x'), (15) 

where the correlator A can be taken to be the usual Dirac delta function. For this prescription 



of quenched noise, however, an analytical approach appears hopeless |32|, |33j. We therefore 
propose to relax this definition and represent the quenched randomness of the soil by a static 
noise rj(x), with correlations 

<r ? (x)7?(x')) = 2D6^(x - x'). (16) 

This form of "columnar" noise, despite being a rather crude approximation, has been previously 
proposed to model soil heterogeneity in cellular automata models of fluvial networks ||34f |. More- 
over, we have found it to be useful for obtaining realistic river networks in numerical simulations 



[35]. In this paper we primarily consider the case of static noise (|lq), corresponding to the limit 
in which the external forcing is constant and the dominant source of noise is the inhomogeneous 
composition of the soil. Results for the case of thermal noise are described in Ref. B. 



3 Statistical solutions 
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3.1 Linear regime 

Even in the absence on any nonlinearity, fundamental conclusions may be drawn from (|l2|). By 
setting g = (i.e., by considering T(h) = \o = const.), we obtain the linear equation 

^ = u\$h + v ± Vlh + 7i, (17) 

which is an anisotropic counterpart of the Edwards- Wilkinson equation |22|, |23| . Let us consider 
for the moment the case of a thermal source of noise r\ with correlations (|l4|). In this case, it can 
be easily shown that the amplitude of the correlation functions along the main directions en and 
ej_ are inversely proportional to the square root of the diffusivities u\\ and v±, respectively. This 
inverse proportionality is well known in the isotropic case (see Ref. [36|, Eqs.(2.19), (2.23), and 



(2.26)). In the anisotropic case, we need only realize that the computation of, say, Cm follows 

— 1/2 

from the correlations of h computed at fixed values of xj_. We then obtain Cm ~ v» and an 



(18) 



analogous result for C±. Thus, in the linear regime (|17|), the ratio of the correlations in the two 
principal directions scales like 

C,| ~ U±, 



In other words, since the preferred direction gives u\\ > u±_ (since the relaxation of material is 
expected to be faster in the direction x\\), the topography is quantitatively rougher, at all scales 
and by the same factor, in the perpendicular direction than in the parallel direction. 



In the case of static noise with correlations (16), it can be shown by dimensional analysis that 



the correlation functions are inversely proportional to the diffusivities. In other words, C ~ v 1 



and Eq. (i8|) is correspondingly changed. The qualitative prediction C± > Cy, however, still 
holds. 

3.2 Nonlinear regime 

In this section we summarize the application of the dynamic renormalization group (DRG) 



|29| , 37, 38, 3S] to our nonlinear model, Eq. fllq). Further details may be found in the Appendix. 

The DRG proceeds in Fourier space by integrating out the fast modes, corresponding to 
large momentum k, over an outer shell A6 _1 < k < A, where A is the upper cutoff in wave 
vector space, and 6 is a rescaling factor. In order to bring the system back to its original size, 
a rescaling is afterwards performed, through the homogeneous transformation 

/i(x||,x ±1 t) =& a - L fc(6" Cx a;||,&- 1 xx,&-* J -t). (19) 

The anisotropy is given by the exponent £x = Cy 1 [compare with Eqs. (|5|) and (||)], and the 
scaling with respect to time, discussed below, is given by the dynamic critical exponent z±. 

After performing this transformation, we are left with an equation with the same form as the 
original one, but with different — renormalized — parameters. The transformation law of these 
parameters under an infinitesimal rescaling b = e di , d£ — > constitutes the flow equations of 
the RG. We are interested in the stable fixed points of these equations, corresponding to scale 
invariant phases in the hydrodynamical limit. 

To first order in the coupling constant A, the flow equations read 

^ = v {] (z ± -2U + \), (20) 
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^QT = "±(^-2), (21) 
^ = A^ ± + 2« ± -2a-^, (22) 

^ = £>(2z]_-2aL-C±-d+l). (23) 
Here we have defined the effective coupling constant 

A= /2 A d - 4 ir d _i, (24) 

with = Sd/{2ir) d and the surface area of a d-dimensional unit hypersphere. The flow 
equations for v± and D are exact to all orders in the perturbation expansion. In the case of 
D, this can be proven to be true for any stochastic equation with a conserved current and a 
non-conserved noise, independently of the details of the current [|40f| . For this results from 
the fact that the nonlinearity is proportional to the external momentum k\\ , and cannot therefore 
renormalize a parameter proportional to k± (see Eq. (|l^)) p9|| . Thus Eqs. (f2l"l) and (p3|) provide 
us with the exact result z± = 2, and with the exact scaling relation 2z± — 2a ± — £j_ = d — 1. 
The dynamic critical exponent z± measures the saturation length of correlations as a function 
of time |ll], 20 1 . However, since the time scales for geomorphologic evolution are many orders of 



magnitude larger than those available for observation, the actual value of z±, appears, at least 
at this point, to be purely a matter of speculation. 

Given (p0|)-(p4|), the effective coupling A flows under rescaling as 

^ = A>-3A), (25) 



where e = 4 — d . The stable fixed points of ( B5h are A| = for d > 4 and A2 = e/3 for d < 4. 
For d > 4 the critical exponents attain their mean-field values a\ = 0, C± = lj an d z^l = 2- On 
the other hand, for d < 4, the critical exponents computed at first order in e are 

a ± = %> C± = 1 + I' Z± = 2 ' (26) 

The physically relevant dimension for erosion in the real world is d = 2. Even though the 
result (^) represents only the first terms in a expansion in powers of e, and is therefore valid only 
approximately for small values of that variable, we can still gain some information by setting 
s = 2. In this case, we obtain the scaling exponents 

5 1 4 a 1 5 . 
a± = - ~ 0.83, C± = -r = -, an = -± = - ~ 0.63. 27 

6 ' S C|| 3 11 Ci 8 V ' 



4 Comparison with field data 

The values (|27| ) predicted for a± and aii are in reasonable agreement with previous measures 
made at small length scales [41-48], where the roughness exponent was reported to be between 
0.70 and 0.85. This observation lends support to the hypothesis of two "universality classes" 
for geomorphological evolution. The first of these classes encompasses topographies evolving 
isotropically at large length scales, and possibly dominated by non-erosive mechanisms such as 
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Figure 2: a) Digital elevation map of an area of the Appalachian Plateau, in Northwest Pennsyl- 
vania. Elevations are given in meters. The spatial resolution is 90 m. b) Averaged height-height 
correlation function C(r) for the landscape in Figure [2]a, where r is oriented in the vertical 
direction of (a). A plot of similar shape, but with smaller values of C(r), is obtained in the 
horizontal case. Logarithms are computed from quantities measured in units of meters. 



internal tectonic stresses. A characteristic of this class is a small roughness exponent which 
is compatible with estimates made from the KPZ equation @. Thus, following the proposal 
of Sornette and Zhang, we can identify the KPZ equation as a description of some universal 
features of the large scale dynamics. On the other hand, for small length scales we expect 
anisotropic effects to be dominant. The anisotropy, which is induced by small-scale inclinations 
of the landscape, would lead to a purely erosive dynamics, which in turn should yield the large 
roughness exponents predicted by our theory. 

If these two universality classes do indeed exist at different length scales, then one should be 
able to find evidence of a crossover from one regime to another in the same piece of topography. 
Specifically, one expects that the correlation function should change from a high a regime 
at short length scales to a low a regime at large length scales. This sort of crossover has 
indeed been reported several times in the literature [41, [4^, |4^, [44|, 45 1. Indeed, Ref. E3] has 
already suggested that the crossover length separates a small-scale erosive regime from large- 



scale tectonic deformation. Ref. [45|, on the other hand, suggests that the crossover separates 
length scales which have had sufficient time to fully develop in a KPZ-like way and those which 
have not. Although it is beyond the scope of this paper to make a definitive argument in favor 
of either of these interpretations, in Figure § we present measurements of our own, from the 
Appalachian Plateau in NW Pennsylvania. Figure || also suggests the presence of a crossover, 
here located at a characteristic scale of about 1 km. We note that in topography depicted in 
Figure ^a, the principal features are deeply eroded channels with a characteristic width of order 
1 km. The long wavelength features, on the other hand, have resulted from tectonic stresses 
associated with the formation of the Appalachian Mountains. Thus, based on the evidence of 
this example, we prefer the interpretation of Ref. [42|. 

In Figure 2 there is no obvious preferred direction, and all of the measurements reported 
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Figure 3: a) Digital map of a submarine canyon off the coast of Oregon, located at coordinates 
44°40' N, 125°45' W. The vertical axis represents the depth z below sea level. The spatial 
resolution is 50 m. b) Mean average profile of the canyon, along the direction x\\. All units are 
given in meters. 



in the literature were either averaged over all directions or the direction of the measurements 
was not reported. Thus, to check the full validity of our results with natural topography that 
has an unambiguous preferred direction, we have analyzed digital bathymetric maps of the 
continental slope off the coast of Oregon. In this case the slope results from the relatively 
abrupt increase in the depth of the seafloor as the continental shelf gives way to the deeper 
continental rise. Figure |3|a shows one portion of this region. Here the main feature of the 
topography is a deep incision called a submarine canyon. In this region, submarine canyons are 
thought to have resulted from seepage- induced slope failure [46|, which occurs when excess pore 



pressure within the material overcomes the gravitational and friction forces on the surface of 
the material, causing the slope to become unstable. Slope instabilities then create submarine 
avalanches, which themselves can erode the slope as they slide downwards. 

In order to make comparisons with our predictions (p7|), we have computed the correlation 
functions C\\ and C±, corresponding, respectively, to the parallel and perpendicular directions of 
the seafloor topography in Fig. ||a. The computation of Cj_ follows naturally from its definition 
but the computation of Cm deserves some comment. The fluctuations measured by Cm must 
be defined with respect to an appropriate average profile. One expects that geologic processes 
other than erosion (e.g., tectonic stresses) are responsible for long- wavelength deformations 
in the parallel direction. Assuming that these deformations are on average constant in the 
perpendicular direction, we may estimate such systematic corrections by computing the mean 
profile along the parallel direction, 



/■/ av (.ci ) = — J dxxh(xpx±), (28) 



where L± is the length of the system in the perpendicular direction. We have plotted h av in 
Fig. ||3. We use it to detrend h by computing the correlation function Cm from the fluctuations 
of the detrended surface h = h — h av (x»). 

Figure |] shows the plots of Cm and Cj_, corresponding to the topography in Fig. ||]a. One sees 
that the least-squares estimates of the roughness exponents, aii ^ 0.67 and a± ~ 0.78, exhibit 
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Figure 4: Height-height correlation functions computed along the parallel (Cm) and perpendicular 
(Cj_) directions for the landscape shown in Fig. ||a. Solid lines are least-squares fits to the scaling 
region. The logarithms are computed from quantities measured in meters. 

a surprisingly good fit to our theoretical predictions (p7|). 

We have also measured Cm and Cj_ in some desert environments. One such example is shown 
in Figure ||, corresponding to an area near Marble Canyon, in Northeast Arizona. In these 
cases we do not obtain conclusive power law scaling, but we almost always find C±/C\\ > 1, as 
predicted by the linear theory. (Indeed, at small scales in this case, we find C± ~ I.8C11) Thus, 
while the example of Figure |3] may be in some sense specialized, one of our main predictions — 
that the topography in the perpendicular direction is rougher than the topography in the parallel 
direction — seems to be of fairly general validity. 



5 Size distribution of avalanches 

Real sloping topography can erode episodically in a series of infrequent events. This episodic 
erosion amounts to a series of avalanches. 

For the case of sloping submarine topography such as that shown in Figure ||a, such avalanches 
can create gravity-driven flows |l4|] of suspended sediment. When these flows finally come to rest, 
the sediment settles out. Then, over geologic time, the sediment solidifies to form sedimentary 



rocks known as turbidites [15]. Partly because these sedimentation events are widely spaced in 
time, individual layers of rock may be associated with each avalanche-like sedimentation event. 
The thickness of these layers may be assumed to be related to the size, or volume of sediment, 
associated with the avalanche that created them. 

Recent empirical studies of turbidite deposits show that in some instances a power-law dis- 
tribution of thicknesses may be observed that extends over nearly two orders of magnitude in 



thickness [16, 17, 18]. In other words, the measurements indicate that the probability P^(A) of 



an avalanche resulting in a deposit of thickness A scales like 

P A (A) ~ A -7 , (29) 
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Figure 5: a) Digital elevation map of an area near Marble Canyon, in Northeast Arizona. 
Elevations are given in meters, and the spatial resolution is 90 m. b) Height-height correlation 
functions computed along the parallel (Cy) and perpendicular (CjJ directions for the landscape 
shown in Fig. ||a. Logarithms are computed from quantities measured in meters. 



where 7 is a characteristic exponent .0 In the best documented cases, 7 is between about 2 and 

2.4 g§ 

Refs. [|17|, |18| suggested that the power-law distribution (|29|) could be the result of a nat- 
ural manifestation of self- organized criticality (SOC) Jl3| , [47| ]. Systems exhibiting SOC, and in 
particular, certain models of sandpiles, exhibit a dynamics dominated by avalanche events, in 
which the number of avalanches of size s scales like the power law 

P s (s) ~ s- T , (30) 

with a characteristic exponent r. Assuming that turbidites result from a series of slumps that 



may be in some way related to the SOC sandpile models, then, as indicated in Refs. |16|, |17], |lg| , 
their size distribution may also be given according to a expression similar to (pCf), One might 
further expect that this distribution could be related to geometric aspects of the surface from 
which the avalanches fall. Our objective in this section, then, is to relate the scaling properties 
of the topography of a sloping surface to the scaling properties of the avalanche size distribution. 

To do so, we follow Refs. pSj , 29] and view our model (|l3|) as a transport equation that 
describes a driven diffusive system, i.e., a sandpile. Under this assumption, we may compute 
the probability distribution of the size of the surface area that relaxes as a result of the interplay 
between the driving force — noise — and the diffusive damping. We assume that this distribution 
is a power law, with a cutoff due to finite-size effects. Moreover, we expect these relaxing surface 
patches to exhibit the same anisotropy as the underlying topography. In other words, as shown 
in Figure [f| we assume that the avalanches result from unstable patches of unit thickness with 

2 Note that our notation differs from Refs. [ jL7| where the cumulative distribution of layers was studied. 
Specifically, their exponent B is equal to, in our notation, 7 — 1. 
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Figure 6: Scaling of an "avalanche patch" over an anisotropic landscape. 



extension proportional to i\\ and £± in the parallel and perpendicular directions, respectively. 
From equations (|5|), (||), and (Q), the self-affine nature of the topography gives £± ~ Thus 

the size s of an unstable patch of extent £y in the parallel direction scales like ^|| + ^" ; and the 

maximum size of an avalanche in a system of parallel extent L scales like L 1+ ^l . The probability 
distribution of avalanche sizes s in a system of finite size L may then be expressed as [^7]] 



P(s,L) = s- r f(-^j. (31) 
Here f(x) is a scaling function such that / is constant for small x < 1 and zero for large x > 1. 



Using equation (|3T[), we can relate r to Oi by means of a scaling argument |34], [48|, |4£|. The 
average size (s) of an avalanche is defined by 



(s) =/" W/ bWj ds - (32) 

Performing the change of variables s = £L 1+ ^l then yields 

( s ) = L^-^+W /y- T /(6d£ (33) 

~ L (2-r)(l+C||). (34) 

On the other hand, since each avalanche is a self-affine patch, as L — > oo the average patch 
should become increasingly elongated in the parallel direction (since Oi < 1). Thus one expects 
(s) ~ L for large L. This relation, together with fl34]), provides us with the result 

(35) 



i + C||' 

which relates the avalanche-size exponent r to the anisotropy exponent £y . 

As noted in Ref. [18, 50 1, we must also address the relationship between the avalanche size 
(or mass) s and the thickness A of the deposited layer. One way to do this is by introducing 
the spreading exponent x such that 

A ~ s x . (36) 



13 



For perfect spreading (i.e., all layers have the same thickness), x = 0, whereas for no spreading 
at all (i.e., all sedimentation events cover the same area A = s/A = const.), x = 1- Empirical 
studies of rock slides, for example, indicate that x ~ 1/3 |5l]], corresponding to self-similar areal 
spreading. From the relation P s (s) = -PA([A(s)])(dA/ds), we find, using equations (31) and 
!§), that 

7 = 1 + —- (37) 

X 

Then from Eq. (|35| ) we obtain 

7 = 1 + -(VM- 08) 



Equation (38) relates the exponents describing surface anisotropy, avalanche size, and spread- 
ing. Using our DRG estimate for the anisotropy exponent, <ji =3/4, and assuming a spreading 
exponent x = 1/3, we find that 

16 , . 

7 = y « 2.3, (39) 



which is in reasonable agreement with the best documented results of Refs. fig , 17, |lq| . On 
the other hand, we may consider equation ( |38|) to be a prediction of the spreading exponent % 
when 7 is measured in a single location but x is unavailable. This could be useful in geological 
applications where one wishes to know the spatial extent of a sequence of turbidite deposits. 



6 Conclusions 

In concluding, it is worthwhile to reflect on the main elements of our theory. Lacking any 
fundamental "equations of motion" for erosion, we have elected to proceed from conservation 
laws and symmetry principles. Thus the principal ingredients of our model are the conservation 
of the eroding material, the presence of a preferred direction for the transport of it, and ran- 
domness in either the landscape or the forcing. Making just these assumptions, we derived an 
anisotropic stochastic equation from which we have extracted both qualitative and quantitative 
predictions. The main qualitative prediction is that eroded topography is rougher in the direc- 
tion across slopes than it is in the direction down slopes. The main quantitative predictions 
are scaling laws for height-height correlation functions. These require additional assumptions 
or restictions concerning the noise and the relevant degree of nonlinearity. Both the qualitative 
and quantitative predictions appear to be in good agreement with measurements made from real 
landscapes. 

We have also included an interpretation of our theoretical model as a driven diffusion equa- 
tion. In this case we have been able to relate the distribution of the sizes of erosion events, or 
"avalanches," to the self-affine scaling that we have predicted for the nonlinear regime of our 
model. The testing of this prediction is beyond the scope of this paper, however, as it would 
require either unusually extensive geological data or an innovative laboratory experiment. 

Our results apply, in principle, to any erosive process with the appropriate lack of symmetry. 
In the usual geological setting, however, the anisotropy applies specifically to a surface of fixed 
inclination which, in turn, implies that our theory should only apply locally, to the relatively 
small scales where the preferred direction of transport is approximately constant. Because the 
anisotropy should vanish at large length scales, we argue that large-scale features of topography 
should be presumably described by a different, isotropic theory, such as the KPZ equation 
@ [11, p^| . We provide evidence, and cite additional results from the literature, that such a 
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crossover indeed exists. We suggest, in line with others |}42fl , that the crossover length separates 
small-scale, externally induced, erosive features of landscapes from large-scale deformations (such 
as those induced by tectonic stresses) of internal origin. 

Finally, we wish to note that, while there is much evidence that landscapes can be self-affine, 
this evidence is rarely unambiguous and certainly not ubiquitous. Our qualitative prediction that 
C± > Cm, on the other hand, appears more robust than any predictions of scaling exponents, or 
even scaling itself. Our results suggest that the coupling of anisotropy to topographic orientation 
may be a fundamental physical property of eroding landscapes. Precisely which length scales 
are relevant to this coupling, however, remains an open question. 

Acknowledgments 

We thank B. Tadic for fruitful discussions and suggestions. R.P.S. acknowledges financial support 
from the Ministerio de Education y Cultura (Spain). The work of D.H.R. was partially supported 
by NSF grant EAR-9706220. 



15 



Appendix 

In this Appendix we develop further details of the renormalization group analysis of Eq. (|l3|), 
where rj is a static noise term, Gaussian distributed, with zero average and correlations 

(^(x^x')) = 2DS^(k - x'). (40) 

Here D is a parameter gauging the strength of the noise. 

In order to proceed, we first Fourier transform the function h, defining 

/i(x,t) = / h{k,uj)e l{k - x -^\ (41) 

Jk 

where we have defined 

d d k r°° duj 



|k|<A (2ir) d J-OC 2tt ' 

The integrals over k are restricted to the upper cutoff A, which plays the role of a lattice spacing, 
or minimum distance in real space. In momentum space, and after performing a few algebraic 
manipulations, Eq. (|l3|) reads 



h(k,u) = G (k,cu)r ? (k,tj) - -kf,G (k,u)) I I fc(q, Si)h(q', Qf)h(k - q - q', w - O - fi'), (42) 



3 ■ JqJq' 

where r/(k, uj) is the Fourier transform of the static noise 7y(x), with correlations 

(T/(k,w)77(q,n)) = 2L>(2^) d+2 (5 (d) (k + q)5(a;)(5(0), 
and the free propagator Go(k, cj) has the form 

Go(k '" > = ^ + ^-^ ' <43) 



In Figure y we have represented Eq. ( |42| ) in terms of Feynman diagrams [^] . 

As we mentioned above, the RG procedure consists of the elimination of the fast modes 
(large wave vectors k), followed by a rescaling of the system back to its original size by means 
of the transformation 

/i(x||,xj_,t) = b a ±h{b-^x\\,b- l x ± ,b~ z H). 

The relevance of the different parameters of the problem (i/m, z^j_, A, and D) can be estimated 
by dimensional analysis. One can check that, under the aforementioned transformation, these 
parameters rescale as 



The coupling constant A will be irrelevant^] when its scaling exponent, A' = b Vx X, is negative, 
where y\ = 2aj_ — 2C,±_ + z± . Selecting the values of z± , a± , and £_L so that the transformations 
of vn, v±, and D are invariant, we obtain 



C° = 1, and a° ± = (4 - d)/2, (45) 



3 In the renormalization group framework, a parameter is irrelevant if it flows to zero, and can therefore be 
neglected in the calculations. 
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Figure 7: a) Diagrammatic expansion of Eq. 
various terms. 



]) in Feynman diagrams, b) Definition of the 



and from them, y\ = 4 — d. The nonlinearity is irrelevant (that is, y\ < 0) when d > d c = 4. This 
defines the critical dimension d c of our system. Above the critical dimension, the nonlinearity 
is negligible, and we recover the scaling exponents (|45|), which are simply given by dimensional 
analysis. On the other hand, below d c the nonlinearity prevails, and we expect fluctuations to 
be dominant and produce nontrivial scaling exponents. 

The RG program is carried out with the help of diagrammatic techniques [ 20 , [37], 36 ] . The 
idea is to graphically iterate (|42| ) to the desired order in A and express the equation in terms 
of effective parameters, which are given as integrals of powers Gq. Afterwards, integration 
of fast modes, by averaging over the noise in the outer shell, and rescaling provides us with 
the renormalized parameters and the consequent flow equations, (|20|)-(|23|), in the limit of an 
infinitesimal transformation. 
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